function [z1,z2,z3,z4] = lpcvd_h2mode_2D_sincos(h,MaxMode)
assert(size(h,1)==size(h,2));
L = size(h,1);
z1 = zeros(MaxMode+1,MaxMode+1);
z2 = z1;
z3 = z1;
z4 = z1;

x_max = pi;
dx = x_max/L;
% dy = L/size(h,2);
x = (0:dx:x_max-dx)';
y = (0:dx:x_max-dx)';

for m = 0:MaxMode
    for n = 0:MaxMode
        phi1mn = fun_phi1mn(m,n,x,y,L);
        phi2mn = fun_phi2mn(m,n,x,y,L);
        phi3mn = fun_phi3mn(m,n,x,y,L);
        phi4mn = fun_phi4mn(m,n,x,y,L);
        z1(m+1,n+1) = sum(sum(h.*phi1mn));
        z2(m+1,n+1) = sum(sum(h.*phi2mn));
        z3(m+1,n+1) = sum(sum(h.*phi3mn));
        z4(m+1,n+1) = sum(sum(h.*phi4mn));
    end
end
    function phi1 = fun_phi1mn(m,n,x,y,L)
        phi1 = (2/L)*sin(2*m*x)*(sin(2*n*y))';
    end

    function phi2 = fun_phi2mn(m,n,x,y,L)
        if (m== 0 || m == L/2) && (n==0||n==L/2)
            coeff = 1/L;
        elseif m~=0 && m~=L/2 && n~=0 && n~=L/2
            coeff = 2/L;
        else
            coeff = sqrt(2)/L;
        end
        phi2 = coeff*cos(2*m*x)*(cos(2*n*y))';
    end

    function phi3 = fun_phi3mn(m,n,x,y,L)
        if n == 0 || n == L/2
            coeff = sqrt(2)/L;
        else
            coeff = 2/L;
        end
        phi3 = coeff*sin(2*m*x)*(cos(2*n*y))';
    end

    function phi4 = fun_phi4mn(m,n,x,y,L)
        if m ==0 || m == L/2
            coeff = sqrt(2)/L;
        else
            coeff = 2/L;
        end
        phi4 = coeff*cos(2*m*x)*(sin(2*n*y))';
    end

end